Self-Organized (Quasi-) Criticality: the Extremal Feder and Feder Model 
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A simple random-neighbor SOC model that combines properties of the Bak-Sneppen and the 
relaxation oscillators (slip-stick) models is introduced. The analysis in terms of branching processes 
is transparent and gives insight about the development of large but finite mean avalanche sizes in 
dissipative models. In the thermodynamic limit, the distribution of states has a simple analytical 
form and the mean avalanche size, as a function of the coupling parameter strength, is exactly 
calculable. 
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I. INTRODUCTION 

Although self-organized criticality (SOC) is perhaps 
not the whole story for the emergence of scale invariance 
in natural systems, certainly it is a promising concept in 
the understanding of how Nature works . The fact that 
still there is not a clear picture of the necessary and/or 
sufRcient mechanisms to create the self-organized critical 
state is reflected in the doubts about whether locally dis- 
sipative systems really presents SOC or have only an ex- 
ponential divergence of the mean avalanche size (s) when 
approaching the conservative limit. The recent result by 
Chabanol and Hakin |^ and by Brock and Grassberger 
1^ stating that the random-neighbor OFC model is not 
critical in the dissipative regime, contradicting a previ- 
ous claim of Lise and Jansen Q , is a clear example of the 
difficulty of making that distinction by using only sim- 
ulation evidence. An instance where this also happened 
is the forest-fire model Q, which has only a very large 
mean avalanche size, not an infinite one like in critical 
systems. 

In this paper, a model is proposed which is similar but 
simpler than the random-neighbor slip-stick model stud- 
ied in H,^. In this model the distribution of states p{E) 
and the mean avalanche size (s), as functions of a dis- 
sipative parameter a, have exactly calculable analytical 
forms (in the infinite system size limit). The analysis in 
terms of branching processes is transparent and gives a 
clear mechanism for the emergence of large but finite (s) 
in systems with local dissipation. The concept of quasi- 
criticality is also discussed as a convenient distinction to 
be made in the classification of the various systems stud- 
ied in the SOC literature. 



II. EXTREMAL FEDER AND FEDER MODEL 

A. The model 

The model is a random-neighbor version of the Feder 
and Feder model |^,^ where it is used an extremal dy- 
namics similar to the Bak-Sneppen model [Q instead of 



the standard global driving. All sites j = 1, . . . ,N have 
a continuous state variable Ej e TZ. At each time step 
the site with maximal value 'fires', resetting its value to 
zero plus a noise term 77. Then, k random 'neighbors' 
(rn) of the firing site have their states incremented by a 
constant a plus a noise term. The choice of neighbors is 
done at the firing instant: the randomness is annealed. 
So, denoting the site with maximal value at instant t as 
E* = ma.x{Ej}, the update rules are: 



EUt- 



Ern{t) 



(1) 



V: 



with 77 being a random variable uniformly distributed in 
the interval [0, e] (the range of e will be discussed later). 
Note that each site receives a different quantity ry. 

Consider the instantaneous density of states p{E, t). It 
is clear that for any E outside the intervals /„ = [{n — 
l)a, (n — l)a -|- ne],n = 1,2, . . ., this density decays to 
zero for long times. These intervals effectively discretize 
the phase space, so it is useful to define the following 
quantities, 



P„. 



{n — l)a-\-ne 



p{E)dE , 



(2) 



(n— 1)q 



max SO that the in- 
max will be obtained 



with rt = 1, 2, . . . , njnax, and e < a/n. 
tervals do not overlap (the integer n. 
latter). The process can be thought of as a transference 
of sites between the intervals J„. At each time step, one 
site is transferred to the interval Ii and, with probability 
kPi, one site is removed from this interval. The aver- 
age flux to the intervals /„ with n > 1 corresponds to 
the probability kPn-i that a neighbor is chosen in the 
previous In-i interval minus the probability fcP„ that a 
neighbor is chosen in the interval /„. The average num- 
ber of sites in each interval is Nn{t) = NPn{t). For long 
times, that is, when the density of states outside the /„ 
intervals goes to zero, one can write 



P^[t + l)=P,{t) + ^[l-kP,{t)] , 

P„(i + 1) = P„(t) + ^ [fcP„-i(t) - kPn{t)] . 



(3) 
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Here, each time step is equal to the update of the maxi- 
mal site and k random neighbors. 

The condition for steady states, P„(t + 1) = Pn{t) = 
P*, gives 



Pi* = l/k 
P* = P 



n-l ; 



(4) 



that is, P* — l/k for all n. But since p{E) is normalized, 
only Umax intervals with P„ of 0(1) can exist. That is. 



n=l 



1 



(5) 



This means that p{E) is composed of k bumps {n — 
1, . . . , Umax = k) and the previous condition for produc- 
ing non-overlapping intervals P„ reads e < a/k. There is 
also a bump of 0(log N/N) (by analogy with the results 
from 1^) situated at the interval Ik+i — [ka, ka + {k + 
l)e]. The other intervals n > k + 1 have P„ of yet smaller 
order (see Fig p. 



B. Avalanches 

An avalanche will be defined as the number of firing 
sites until an extremal site value falls bellow the thresh- 
old Eth — 1. Note that the first site of an avalanche 
(the 'seed') always has E < 1 but it counts as a firing 
site. So, if a seed produces no supra-threshold sites ('de- 
scendants'), this counts as an avalanche of size one. This 
definition of avalanches agree with that used in the stud- 
ies of relaxation oscillator models. 

In these random neighbor models, an avalanche can 
be identified as a branching process where an active site 
produces k new sites, each one having a probability p of 
being active (a 'branch') and a probability 1 — p of being 
inactive (a 'leave'). The branching ratio a — kp measures 
the probability that a firing site produces another firing 
site. Now, consider an avalanche which has terminated 
after s sites have fired. This avalanche is composed by 
one seed and s — 1 descendants. But the average number 
of descendants produced by s firing sites is as. Then, in 
average, the relation 



{s)-l = a{s) 



must hold, which leads to 



is) 



l-a 



(6) 



(7) 



This result also can be obtained from percolation theory 
in the Bethe lattice. Note that a = a{t — > oo) refers to 
the stationary value of the branching ratio: during the 
transient, cr(t) is a function of i. Indeed, the evolution of 
a toward cr(cxD) corresponds to the 'self-organization' of 
the system. 



C. The Q = 1/fc case 

In the case a = 1/fc the calculation of cr is trivial. The 
fc-th bump, which starts at (fc — l)a, must lie bellow the 
threshold Eth = 1 (if not, the system is supercritical). 
Then, e must satisfy the condition (fc — l)/fc + ke < 1, 
that is. 



e < l/fc2 



(8) 



For the standard k — 4 neighbor case this reads e < 
0.0625. This condition also implies that neighbors per- 
taining to the other bumps do not contribute to cr, that is, 
cannot fire when receiving a maximal contribution a + e. 
Now, since all the neighbors pertaining to the fc-th bump 
receive at least the quantity a — 1/fc, they are determin- 
istically transformed into active sites. Thus, the average 
number of descendants of a firing site is 



a = kP^ = k X y = 1 , 



(9) 



This corresponds to a critical branching process. It is 
know that in this case the system presents an infinite (s) 
(see Eq. (|^)) and a pure power law P(s) — cs^^/"^ for the 
distribution of avalanche sizes ||] . 



D. Results for general a 

For a < 1/fc, the distribution of states p{E) must be 
known. But it is clear that if ka = 1 — S then inevitably 
cr < 1 (even for very small S > 0), since some sites per- 
taining to the fc-th bump may not receive a sufficient 
contribution to make them active (see Eq. ( p^ below). 
Thus, any value a < ac — 1 / k is subcritical. 

The calculation oip{E) is very simple. In the station- 
ary state, a site pertaining to the n-th bump has state 
E = (n — l)a -\- Zm where z„ is the sum of n random 
variables uniformly distributed in the interval [0, e]. The 
distribution p{zn) may be calculated from 



Pizi) 

PiZn+l) 



ie(zi)e(e-zi) 



(10) 



dZndzi p{Zn)p{zi)S{Zn + Zi - Zn+l) 



For the fc = 4 case, 

p(z2) = e-2 [Z2e(z2)e(e-Z2) 

+ (2e-Z2)e(z2 -e)e(2e-Z2)] 



(11) 



P{z3) = 



^e(z3)e(e-z3) 



3ez3 



3ez3 



9e^ 



e(z3-e)e(2e-Z3) 



e(z3 - 2e)e(3e - Z3) 



2 



p{z4) = 



r^3 



6 



■e(z4)e(e-z4) 



^4 , 0,^2 0,2 2e' 



2 , 2ezi - 2e'z4 + — ) 6(24 - e)e(2e - Z4) 
^ + 2ex^ - 2e^x + ^ ) 6(^4 - 2e)e(3e - Z4,) 



+ — e(z4-3e)e(4e-Z4) 
6 



with the shorthand x = (4e — Z4). Each bump with label 
n in the distribution p{E) starts at En — {n — l)a, being 
proportional to p{zn) (the constant of proportionality is 
just l/k). In Fig. |l|, the distribution p{E) is compared 
with simulation results for a system with = lO** sites, 
a = 0.235 and e = 0.05. 



if they receive a quantity 01+77 with r] > S—z. This occurs 
with probability P(77 > S — z) = l — {S — z)/€. Then, these 
sites contribute to the branching ratio with the quantity 



piz) 1 - 



S^z 



dz . 



(13) 



The total branching ratio is then 



a = a' + a" = 1 



p(z) dz 



61 1 C 

- I p{z) dz — / zp{z)dz 
e Js-e e Js-e 



(14) 



where it was used the fact that p{z)dz — 1. 

Since p{z) has a simple piece-wise polynomial form (see 
Eq. (pT|)) the calculation of (s) is straightforward and the 
result is presented in Fig. || along with simulation results 
for the = 4 case with e = 0.05. 
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FIG. 1. Distribution of states p{E) for fe = 4 and e = 0.05: 
theoretical (solid) and simulation (circles) with A'^ = 10* sites. 
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The branching ratio a is calculated as follows. All the 
sites that can be activated are in the fc-th bump. When 
they are hit, sites with E > \ — a activates deterministi- 
cally. In terms of the re-scaled variable Zk = E—{k — l)a, 
this condition refers to sites with Zk > 5 = 1 — ka. They 
contribute to the branching ratio with the quantity ct'. 



a' = k 



p{E)dE 



S+a 



p{z) dz , 



(12) 



where z = z^. 

Sites with z < 5 — e cannot be activated and does not 
contribute to a. Sites with 5 — t < z < 5 can be activated 



FIG. 2. Mean avalanche size (s) as a function of parameter 
a. Theoretical (solid) and simulations with A' — 10* sites for 
noisy EFF model (circles) with e = 0.05 and noiseless EFF 
model (triangles) with e = 0.2. These e values are chosen 
such that the last interval (74) has the same length in both 
models. 

For 5 < e, that is, Uc — a < e/fc, the form assumed by 
a is particularly simple, since p{z) = Ce~ 
interval (C is a numerical constant). Then, 



C 



z''-\S-z)dz 



3 



c 



k{k + l) 



k+1 



(15) 



Since S = 1 — ka = k{ac — a), this means that (s) diverges 
as (ttc — a)^'^^^, see Eq. (|^). For example, for k — A and 
€ = 0.05, (s) = 120 ah-eady for a = 0.2375. Curiously, 
this behavior is similar to the (etc— a)^*^ divergence found 
in the standard random neighbor FF model 



E. The 'noiseless' EFF model 

It is instructive to compare this behavior with that of 
a simpler FF model where the firing rule is the same, 
E*{t + 1) — T] G [0, e] but the neighbor update rule 
is noiseless, Emit + 1) = Emit) + a. Then, p{E) as- 
sumes the form of k rectangular bumps with p{zn) — 
e~^9(z„)9(e — z„). In this case, here called 'noiseless' 
EFF model, the branching ratio is ct = for 5 > e and 
a — 1 ~ 6/e for < S < e. In contrast with the previous 
model, large avalanches only occur when a is very close 
to ac (see Fig. |). 



III. COMPARISON WITH THE 
RANDOM-NEIGHBOR BAK-SNEPPEN MODEL 

The Bak-Sneppen model [Q consists of N sites car- 
rying real variables x (fitness) that follow an extremal 
dynamics with the rules: 



X* = minjxj} , 
x*it + l)=7], 

Xrnit + 1) = ?7 , 



(16) 



where rj is an uniformly random variable in the interval 
[0, 1]. To facilitate the comparison with the extremal FF 
model, one makes a change of variables to E = 1 — x, 
the BS neighbor update rule may be written as Em — 
1] G [0,^] (the Bak-Sneppen choice is A = 1). Thus, the 
parameter A is related to the effective coupling strength 
between neighbors. 

In the iV — > cx) limit, its stationary distribution of 
states is p{E) = (fc_\)^ 6( ^'^'^^^^ - x). Note that, here, 
k is not the number of neighbors q but the number of 
tentative branches k = 1 + q since, in the BS model, the 
firing site is also a tentative branch which may contribute 
to a. Thus, the random neighbor BS model presents, in 
the stationary state, a natural threshold at 



Xc = {k-l)A/k. 



(17) 



All firing ('mutating') sites have E above Ac, so Bak- 
Sneppen argued that it is natural to choice Ac as the 
threshold that defines avalanches. With this choice, the 
distribution p{s) is a power law. 

Like the Bak-Sneppen model, a 'natural' threshold 
Ac (a) also appears in the extremal FF system. Suppose 



that the initial condition p{E, t = 0) is uniform in the 
interval [0,1]. Then, define the 'gap' at time t as the 
minimal extremal site chosen up to t. During the tran- 
sient, the gap evolves toward the value 



Ac(a) = (k - l)a + ke 



(18) 



Now, if we define an avalanche as the return time between 
two crosses of Ac by the time series of extremal values, 
like is done by Bak and Sneppen, then the system will 
stay always critical since Ac is a function of a and all is 
the same when a is changed. Indeed, any choice for A in 
the interval [Ac (a), -Be (a) — ka] makes the A-avalanches 
critical. 

Thus, the Bak-Sneppen choice for the threshold A as 
the point where p{E) drops, although natural, is indeed 
a fine tuning operation! It corresponds, in the extremal 
FF model, to choose the proper Ac(a) (or, alternatively, 
Ec{a) — ka) for defining avalanches. Then, in the EFF 
model, instead of fixing Eth = 1 and fine tuning the cou- 
pling parameter to adEth) = 1/fc, the parameter a could 
be fixed and the threshold could be fine tuned to Ec{a) 
(or Ac(q;)). The fine tuning is incorporated (hidden?) in 
the definition of avalanches! 

Now, for comparing between the coupling strength in 
the FF model (in the limit e ^ 0) and the BS model, 
calculate the average new value assumed by a neighbor. 



BS: {E{t + 1)) 
FF: {E{t + 1)) 



A/2, 

{E) + 



ka/2 



(19) 



since {E) = (fc - 1)q;/2 for the FF model. Thus, to 
see the strong analogy between the extremal FF and the 
BS models, it is natural to write A — ka (remember, 
A is an arbitrary parameter which controls the effective 
coupling between sites in the BS model). Then, Eq. (17) 
is Ac = (fc — l)a which has the same form of Eq. (IS) 
(with e ^ 0). This means that there is a well defined 
dependence of Ac on the parameter a which controls the 
coupling between sites: the definition of Ac, although 
natural, is not robust to changes in the coupling strength. 
The issue here is not robustness but automatic fine tuning 
of the threshold Ac (see bellow). 



IV. ON SOC DEFINITIONS 

The idea of self-organized criticality present in the lit- 
erature embodies two distinct properties. The term self- 
organized refers to the fact that there exist a parame- 
ter (cr(t)), which controls the avalanche behavior, whose 
value is not fixed a priori like, for example, in standard 
percolation, but evolves during a transient phase toward 
an stationary value a{oo). Indeed, this time dependence 
should be written as a{p{E,t)), that is, cr is a func- 
tional of the time varying distribution of states p{E, t) 
and a{oo) = a(j){E, oo)). 

The evolution of p{E, t) toward the non-equilibrium 
steady-state p{E, oo) is akin to the transient relaxation 



4 



in equilibrium systems: any initial condition leads to the 
same stationary state, thus, the same value for (t(oo). 
However, this robustness to initial conditions and per- 
turbations ('dynamical stability') should not be mistaken 
as parameter robustness ('structural stability'). Then, a 
second and distinct characteristic recurrently claimed for 
SOC systems is structurally stable criticality |^ ||^: there 
exists a finite parameter range where, after the transient, 
the system is critical (that is, (7{oo) = ctc). For this sit- 
uation, no fine tuning is needed to find a system with 
critical behavior, only some gross tuning. This kind of 
structurally stable criticality will be called type-I SOC. 

Now, it seems that all the SOC models in the literature, 
at least in the random-neighbor versions or other analyt- 
ically tractable situations, are not robust to changes in 
the coupling strength between sites. For example, it is 
well known that the sand pile model is not critical if there 
exists a non vanishing probability of a grain to evaporate 
during an avalanche. This probability plays a similar role 
of the quantity 1 — fca in the present model. 

Structurally stable criticality is depicted in Fig. ||(a). 
In this case, by changing the coupling parameter a, there 
is a finite range where a assumes the critical value Uc = 1. 
Standard criticality, then, may be viewed as the fine 
tuning dependence of cr on a represented in Fig. ^(c). 
However, a third possibility is a strong nonlinear depen- 
dence of cr{a) (Fig ^(b)) like the observed in the noisy 
extremal FF model and in the standard random neigh- 
bor FF and OFC models studied in The system is 
almost critical in a large parameter region: this will be 
called generic quasi- criticality . In standard critical phe- 
nomena parlance, the system has a very large 'critical 
region'. The importance of this characterization is that 
various systems examined in the SOC literature, previ- 
ously beheld as having true generic criticality, are now 
being recognized as having only generic quasi-criticality 
in the sense defined above. 

A candidate for type-I SOC is the two-dimensional 
OFC model Of course, models defined in lattices 

are hard to analyze. It is not clear if the other factors 
present in that models (spatial structure, large probabil- 
ity of re- visitation of sites during an avalanche, boundary 
conditions etc.) create a truly generic critical state for 
a finite interval in the coupling parameter (a) space or 
only a generic quasi-critical behavior. It seems reason- 
able to ask, however, why these arbitrary factors should 
conspire to give exactly a critical state 17(00) = 1 instead 
of only amplifying the tendency toward large avalanches 
already present in the random-neighbor version. If the 
latter view is correct, then a conjecture can be made that 
a necessary condition for lattice models (like 2D OFC) 
to present an apparent generic SOC behavior is that the 
corresponding random-neighbor version presents strong 
quasi-criticality behavior. 

It also has been argued that, by analogy with the ex- 
tremal FF model, the Bak-Sneppen model should not 
be viewed as having structural stability. The choice of 
the threshold Ac for defining avalanches, although natu- 



ral, is similar to the fine dependence Ac(ct) or Ec{a) in 
the present model. Ironically, the random-neighbor BS 
model is even not 'self-organized' in the sense defined 
here: with A = Ac fixed, the branching ratio is always 
a ~ 1 and does not evolve during the transient phase. In 
contrast with the EFF model, the probability of a neigh- 
bor site to fire does not depend on its previous state, so 
that there is no coupling between a and P{E, t). 




a 



FIG. 3. a) Structurally stable criticality: the value of pa- 
rameter a is critical on a finite range of the system parameter 
a. b) Generic quasi-criticality (noisy EFF model): the depen- 
dence of cr{a) is non-linear near ctc c) Standard criticality 
(noiseless EFF model): the system parameter a must be fine 
tuned for obtaining ct ~ ctc. 

Although type-I SOC (or SSC, structurally stable crit- 
icality) seems to be hard to find, another possibility con- 
sists in systems where the coupling strength a is not a 
fixed parameter but turns out a slow dynamical variable. 
That is, although criticality requires a fine tuning of a, 
this is done automatically by some self-tuning mecha- 
nism. This wiU be called type-II SOC or STC (self-tuned 
criticality). For example, suppose that the coupling ay- 
evolves as 



a,jit + l) = il + ni)a,jit) -n2A,{t) , 



(20) 



with fii < ^2, both small quantities, and Ai{t) = 1 if site 
j has fired at time t and Ai{t) = otherwise. Thus, a su- 
percritical state will present a lot of big avalanches which 
shall decrease the average coupling toward the subcritical 
regime. On the other hand, if the system is subcritical, 
the average (a^ ) grows. This dynamics for the couplings 
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has a stationary state P{a) that seems to produce a ro- 
bust critical state . Note that the coupHng dynamics 
is local (depends only on site i) and does not require a 
fine timing in the parameters Qi and ^2- 

The motivation for this dynamics may be (in the con- 
text of neural networks) an activity dependent depression 
in synaptic efhcacy due to synaptic fatigue. In slip-stick 
context, this dynamics may correspond to some plastic 
change in the coupling mechanism between sites (for ex- 
ample, rupture and recover of contacts between neigh- 
bors). Systems with type-II SOC, however, have not yet 
been fully explored in the literature (see |11[). 
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V. CONCLUSIONS 



A class of extremal slip-stick models has been intro- 
duced. The FF model has been studied in the iV — > oo 
limit; the model presents a clear mechanism for produc- 
ing a strong divergence on the average size of avalanches 
like that found in another more complicated models [^,|| . 
As possible extensions, one could examine finite dimen- 
sional versions of these extremal models and the use of 
other update rules (say, OFC rules). Preliminary results 
shows synchronization phenomena for the ID case 

From a practical viewpoint, that is, for explaining 
generic scale invariance in Nature, structurally stable 
criticality or generic quasi-criticality are almost identical 
(up to a cutoff size always present in natural systems). 
It is, however, of interest to determine what are the ba- 
sic ingredients for producing generic quasi-criticality in 
the models examined in the SOC literature. The simple 
mechanism devised in this work suggests that, if true ro- 
bust criticality (type-I SOC) is not easy to be achieved, 
generic quasi-criticality certainly is. Perhaps, true type- 
I SOC do not exist, only generic quasi-criticality and, 
perhaps, type-II (self-tuned) criticality. 

Aknowledgments: The author thanks C. P. C. 
Prado and S. R. A. Salinas for useful suggestions, N. 
Dhar for illuminating discussions about the SOC con- 
cept, J. F. Fontanari and R. Vicente for commenting the 
manuscript and FAPESP for financial support. 



[1] P. Bak, How Nature Works (Copernicus, New York, 
1996). 

[2] M-L. Chabanol and V. Hakin, Phys. Rev. E 56 R2343 
(1997). 

[3] H-M. Broker and P. Grassberger, Phys. Rev. E 56 3944 
(1997). 

[4] S. Lise and H. J. Jensen, Phys. Rev. Lett. 76 2326 (1996). 
[5] G. Grinstein, in Scale Invariance, Interfaces and Critical- 
ity, Vol. 344 of NATO Advanced Study Institute, Series 



6 



